# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools", "sdmTMB", "sdmTMBextra", "terra", "mapplots",
"viridis", "visreg", "modelr", "future", "kableExtra", "ggh4x", "patchwork")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Import some plotting functions
# Source code for map plots
# You need:
# devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")
#remotes::install_github("pbs-assess/sdmTMBextra", dependencies = TRUE)
library(ggsidekick)
theme_set(theme_sleek())
# Set path
home <- here::here()Year-length-effects_Feeding-ratio-total
Background
We want to understand how the effects of the length of the predator (cod) influences the total amount of food it has in its stomach and if this has changed over the course of the time series (1963-2022). Using two different models, we find that this effect depends on model configuration. In the first model (M1f below), the effect of time (Year) has a ….
We test four different ways of modellkng this interaction between Year and predator lenngth (pred_length_sc) and compare
Depending on model, either the predator length effect varying with Year (M1, pred_length_sc) or the Year effect (in M2) is sucking up the trend (pardon my lack of correct lingo).
Load libraries
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/main-analysis/03-fit-diet-models_cache/html"))Read stomach data
df <- read_csv(paste0(home, "/data/clean/stomachs.csv")) |>
mutate(depth_sc = (depth - mean(depth))/sd(depth),
herring_sc = (herring - mean(herring))/sd(herring),
saduria_sc = (saduria - mean(saduria))/sd(saduria),
sprat_sc = (sprat - mean(sprat))/sd(sprat),
other_invert_sc = (other_invert - mean(other_invert))/sd(other_invert),
other_sc = (other - mean(other))/sd(other),
other_fish_sc = (other_fish - mean(other_fish))/sd(other_fish),
benth_fish_sc = (benth_fish - mean(benth_fish))/sd(benth_fish),
Year_f = as.factor(Year),
month_f = as.factor(Month),
ices_rect = as.factor(ices_rect),
pred_length_sc = (pred_length - mean(pred_length)) / sd(pred_length))
glimpse(df)Rows: 126,370
Columns: 38
$ pred_ID <chr> "DK_2007_11_11_160_40G6_13", "DK_2007_11_11_162_40G6_1…
$ herring <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ saduria <dbl> 0.00, 0.00, 0.00, 0.94, 0.00, 22.37, 0.00, 2.93, 1.43,…
$ sprat <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ other <dbl> 36.67, 0.00, 0.00, 0.00, 42.19, 85.10, 0.00, 76.07, 53…
$ other_invert <dbl> 0.00, 0.29, 0.00, 0.03, 0.03, 0.00, 0.00, 0.33, 0.00, …
$ other_fish <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ benth_fish <dbl> 0.00, 0.12, 0.00, 0.00, 0.00, 0.52, 0.00, 0.00, 0.00, …
$ fr_tot <dbl> 2.860374e-02, 9.761905e-03, 0.000000e+00, 5.549199e-04…
$ fr_sad <dbl> 0.0000000000, 0.0000000000, 0.0000000000, 0.0005377574…
$ fr_spr <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ fr_her <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ fr_other <dbl> 2.860374e-02, 0.000000e+00, 0.000000e+00, 0.000000e+00…
$ fr_other_invert <dbl> 0.000000e+00, 6.904762e-03, 0.000000e+00, 1.716247e-05…
$ fr_benth_fish <dbl> 0.0000000000, 0.0028571429, 0.0000000000, 0.0000000000…
$ fr_other_fish <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ Year <dbl> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, …
$ Month <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11…
$ day_of_year <dbl> 315, 315, 315, 315, 315, 315, 315, 315, 315, 315, 315,…
$ pred_weight <dbl> 1282, 42, 38, 1748, 1636, 1988, 1860, 2714, 2270, 52, …
$ pred_length <dbl> 54, 17, 16, 57, 59, 62, 58, 67, 65, 17, 18, 15, 18, 17…
$ lat <dbl> 55.75, 55.75, 55.75, 55.75, 55.75, 55.75, 55.75, 55.75…
$ lon <dbl> 16.5, 16.5, 16.5, 16.5, 16.5, 16.5, 16.5, 16.5, 16.5, …
$ ices_rect <fct> 40G6, 40G6, 40G6, 40G6, 40G6, 40G6, 40G6, 40G6, 40G6, …
$ X <dbl> 594.1508, 594.1508, 594.1508, 594.1508, 594.1508, 594.…
$ Y <dbl> 6179.275, 6179.275, 6179.275, 6179.275, 6179.275, 6179…
$ depth <dbl> 60.06, 60.06, 60.06, 60.06, 60.06, 60.06, 60.06, 60.06…
$ depth_sc <dbl> -0.4524581, -0.4524581, -0.4524581, -0.4524581, -0.452…
$ herring_sc <dbl> -0.1267861, -0.1267861, -0.1267861, -0.1267861, -0.126…
$ saduria_sc <dbl> -0.160259566, -0.160259566, -0.160259566, -0.096552035…
$ sprat_sc <dbl> -0.1886723, -0.1886723, -0.1886723, -0.1886723, -0.188…
$ other_invert_sc <dbl> -0.2556360, -0.2001578, -0.2556360, -0.2498969, -0.249…
$ other_sc <dbl> 3.09833934, -0.06401689, -0.06401689, -0.06401689, 3.5…
$ other_fish_sc <dbl> -0.02739089, -0.02739089, -0.02739089, -0.02739089, -0…
$ benth_fish_sc <dbl> -0.06023454, -0.05066107, -0.06023454, -0.06023454, -0…
$ Year_f <fct> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, …
$ month_f <fct> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11…
$ pred_length_sc <dbl> 1.222301, -1.286201, -1.353999, 1.425693, 1.561288, 1.…Filter
# remove empty stomachs
df_noz <- df |>
filter( fr_tot > 0 )filter: removed 43,832 rows (35%), 82,538 rows remaining# remove na day of Year if that is used as a predictor
df_noz_doy <- df_noz |>
filter( !is.na(day_of_year) ) |>
mutate(doy_sc = (day_of_year - mean(day_of_year))/sd(day_of_year)) # read_csv(paste0(home, "/data/clean/stomachs.csv")) |> summarise(mean = mean(pred_length))filter: removed 422 rows (1%), 82,116 rows remainingmutate: new variable 'doy_sc' (double) with 280 unique values and 0% NA# reduce number of Years to reduce fitting time
df_noz_doy_sel <- df_noz_doy |> filter( Year %in% c(1982:2022))filter: removed 42,762 rows (52%), 39,354 rows remainingmesh_nozdoy_sel <- make_mesh(df_noz_doy_sel, c("X", "Y"), cutoff = 6)
# missing Years
my <- min(df_noz_doy_sel$Year):max(df_noz_doy_sel$Year)
missing_Years <- my[!my %in% unique(df_noz_doy_sel$Year)]
# df_noz |>
# filter(Year %in% c(2013:2016)) |>
# ggplot(aes(pred_length,pred_weight, color = Year)) +
# geom_point()
#
# df_noz |>
# filter(Year %in% c(2013:2016)) |>
# ggplot(aes(pred_length,fr_tot, color = Year)) +
# geom_point()We suspect that the amount of food in relation to cod weight in stomachs is decreasing and depends on predator length:
df_noz_doy_sel |>
group_by(Year) |>
summarise(mean_fr_tot = mean(fr_tot),
sd_fr_tot = sd(fr_tot)) |>
ggplot(aes(Year, mean_fr_tot)) +
geom_line() +
geom_ribbon(aes(ymin = mean_fr_tot - sd_fr_tot, ymax = mean_fr_tot + sd_fr_tot), alpha = 0.3) +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) group_by: one grouping variable (Year)summarise: now 39 rows and 3 columns, ungrouped
df_noz_doy_sel |>
filter(Year %in% seq(1982,2022, by = 10)) |>
ggplot(aes(pred_length, fr_tot, color = Year)) +
geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k=3))filter: removed 35,195 rows (89%), 4,159 rows remainingWarning: The following aesthetics were dropped during statistical transformation: colour
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
Problem based on previous models (M1 and M2)
Year modeled as a fixed effect factor (M1) or a time varying intercept (M2)
M1 <-
sdmTMB(
data = df_noz_doy_sel,
formula = fr_tot ~ 0 + as.factor(Year) + s(doy_sc, bs = "cc") + depth_sc,
time_varying = ~ 0 + pred_length_sc, # implies time varying random walk intercept
extra_time = missing_Years, # to fill in empty Years
time = "Year", # for spatiotemporal and time_varying
mesh = mesh_nozdoy_sel,
spatial = "off",
family = Gamma(link = "log"),
)
M2 <-
sdmTMB(
data = df_noz_doy_sel,
formula = fr_tot ~ 0 + s(doy_sc, bs = "cc") + depth_sc,
time_varying = ~ 1 + pred_length_sc,
extra_time = missing_Years, # to fill in empty Years
time = "Year", # for spatiotemporal and time_varying
mesh = mesh_nozdoy_sel,
spatial = "off",
family = Gamma(link = "log"),
)
sanity(M1)✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largesanity(M2)✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeqqnorm(residuals(M1))
qqline(residuals(M1))
qqnorm(residuals(M2))
qqline(residuals(M2))
M1Spatiotemporal model fit by ML ['sdmTMB']
Formula: fr_tot ~ 0 + as.factor(Year) + s(doy_sc, bs = "cc") + depth_sc
Mesh: mesh_nozdoy_sel (isotropic covariance)
Time column: Year
Data: df_noz_doy_sel
Family: Gamma(link = 'log')
coef.est coef.se
as.factor(Year)1982 -3.66 0.09
as.factor(Year)1983 -3.77 0.07
as.factor(Year)1984 -3.99 0.11
as.factor(Year)1985 -3.58 0.13
as.factor(Year)1986 -3.50 0.11
as.factor(Year)1987 -3.87 0.13
as.factor(Year)1988 -3.76 0.15
as.factor(Year)1989 -3.72 0.13
as.factor(Year)1990 -3.49 0.12
as.factor(Year)1991 -3.75 0.14
as.factor(Year)1992 -3.42 1.25
as.factor(Year)1993 -3.99 0.14
as.factor(Year)1994 -3.82 0.19
as.factor(Year)1995 -4.22 0.17
as.factor(Year)1996 -4.32 0.17
as.factor(Year)1997 -4.13 0.24
as.factor(Year)1998 -3.95 0.21
as.factor(Year)1999 -4.19 0.14
as.factor(Year)2000 -4.38 0.14
as.factor(Year)2001 -4.22 0.17
as.factor(Year)2002 -4.33 0.20
as.factor(Year)2003 -4.18 0.15
as.factor(Year)2004 -4.46 0.14
as.factor(Year)2005 -4.10 0.13
as.factor(Year)2006 -4.42 0.10
as.factor(Year)2007 -4.26 0.09
as.factor(Year)2008 -4.36 0.09
as.factor(Year)2009 -4.40 0.09
as.factor(Year)2010 -4.27 0.14
as.factor(Year)2011 -3.81 1.25
as.factor(Year)2012 -4.51 0.09
as.factor(Year)2013 -4.30 0.05
as.factor(Year)2014 -4.41 0.18
as.factor(Year)2015 -4.18 0.15
as.factor(Year)2016 -4.44 0.11
as.factor(Year)2017 -4.46 0.10
as.factor(Year)2018 -4.27 0.09
as.factor(Year)2019 -4.53 0.16
as.factor(Year)2020 -4.26 0.11
as.factor(Year)2021 -4.38 0.06
as.factor(Year)2022 -4.51 0.12
depth_sc 0.00 0.01
Smooth terms:
Std. Dev.
sds(doy_sc) 0.05
Time-varying parameters:
coef.est coef.se
pred_length_sc-1982 0.36 0.04
pred_length_sc-1983 0.28 0.03
pred_length_sc-1984 0.23 0.04
pred_length_sc-1985 0.09 0.04
pred_length_sc-1986 0.10 0.04
pred_length_sc-1987 0.06 0.04
pred_length_sc-1988 0.16 0.05
pred_length_sc-1989 0.04 0.05
pred_length_sc-1990 0.10 0.04
pred_length_sc-1991 0.26 0.05
pred_length_sc-1992 0.07 0.15
pred_length_sc-1993 -0.12 0.07
pred_length_sc-1994 -0.04 0.08
pred_length_sc-1995 0.00 0.05
pred_length_sc-1996 -0.12 0.07
pred_length_sc-1997 0.26 0.08
pred_length_sc-1998 0.10 0.08
pred_length_sc-1999 -0.01 0.04
pred_length_sc-2000 -0.12 0.06
pred_length_sc-2001 -0.03 0.05
pred_length_sc-2002 0.21 0.08
pred_length_sc-2003 -0.05 0.06
pred_length_sc-2004 0.02 0.04
pred_length_sc-2005 0.06 0.05
pred_length_sc-2006 0.31 0.04
pred_length_sc-2007 0.23 0.03
pred_length_sc-2008 0.43 0.03
pred_length_sc-2009 0.42 0.03
pred_length_sc-2010 0.28 0.04
pred_length_sc-2011 0.39 0.15
pred_length_sc-2012 0.49 0.05
pred_length_sc-2013 0.37 0.03
pred_length_sc-2014 0.64 0.14
pred_length_sc-2015 0.67 0.09
pred_length_sc-2016 0.39 0.06
pred_length_sc-2017 0.12 0.05
pred_length_sc-2018 0.48 0.05
pred_length_sc-2019 0.15 0.07
pred_length_sc-2020 0.46 0.05
pred_length_sc-2021 0.02 0.04
pred_length_sc-2022 0.09 0.06
Dispersion parameter: 0.83
Matérn range: 15.36
Spatiotemporal IID SD: 0.54
ML criterion at convergence: -120570.446
See ?tidy.sdmTMB to extract these values as a data frame.M2Spatiotemporal model fit by ML ['sdmTMB']
Formula: fr_tot ~ 0 + s(doy_sc, bs = "cc") + depth_sc
Mesh: mesh_nozdoy_sel (isotropic covariance)
Time column: Year
Data: df_noz_doy_sel
Family: Gamma(link = 'log')
coef.est coef.se
depth_sc 0 0.01
Smooth terms:
Std. Dev.
sds(doy_sc) 0.05
Time-varying parameters:
coef.est coef.se
(Intercept)-1982 -1.87 0.03
(Intercept)-1983 -1.88 0.03
(Intercept)-1984 -1.89 0.03
(Intercept)-1985 -1.86 0.03
(Intercept)-1986 -1.85 0.03
(Intercept)-1987 -1.87 0.03
(Intercept)-1988 -1.87 0.04
(Intercept)-1989 -1.87 0.04
(Intercept)-1990 -1.87 0.04
(Intercept)-1991 -1.90 0.04
(Intercept)-1992 -1.94 0.04
(Intercept)-1993 -1.98 0.04
(Intercept)-1994 -2.00 0.04
(Intercept)-1995 -2.04 0.04
(Intercept)-1996 -2.07 0.04
(Intercept)-1997 -2.08 0.04
(Intercept)-1998 -2.09 0.04
(Intercept)-1999 -2.11 0.04
(Intercept)-2000 -2.13 0.04
(Intercept)-2001 -2.14 0.04
(Intercept)-2002 -2.15 0.04
(Intercept)-2003 -2.15 0.04
(Intercept)-2004 -2.16 0.04
(Intercept)-2005 -2.15 0.03
(Intercept)-2006 -2.17 0.03
(Intercept)-2007 -2.16 0.03
(Intercept)-2008 -2.18 0.03
(Intercept)-2009 -2.19 0.03
(Intercept)-2010 -2.19 0.03
(Intercept)-2011 -2.19 0.04
(Intercept)-2012 -2.20 0.03
(Intercept)-2013 -2.18 0.02
(Intercept)-2014 -2.18 0.03
(Intercept)-2015 -2.18 0.03
(Intercept)-2016 -2.20 0.03
(Intercept)-2017 -2.20 0.03
(Intercept)-2018 -2.19 0.03
(Intercept)-2019 -2.19 0.03
(Intercept)-2020 -2.19 0.03
(Intercept)-2021 -2.20 0.03
(Intercept)-2022 -2.22 0.04
pred_length_sc-1982 0.35 0.04
pred_length_sc-1983 0.28 0.03
pred_length_sc-1984 0.25 0.04
pred_length_sc-1985 0.10 0.04
pred_length_sc-1986 0.11 0.04
pred_length_sc-1987 0.05 0.04
pred_length_sc-1988 0.16 0.05
pred_length_sc-1989 0.05 0.05
pred_length_sc-1990 0.12 0.04
pred_length_sc-1991 0.26 0.05
pred_length_sc-1992 0.07 0.15
pred_length_sc-1993 -0.13 0.07
pred_length_sc-1994 -0.03 0.08
pred_length_sc-1995 -0.01 0.05
pred_length_sc-1996 -0.14 0.07
pred_length_sc-1997 0.26 0.08
pred_length_sc-1998 0.09 0.08
pred_length_sc-1999 -0.01 0.04
pred_length_sc-2000 -0.12 0.06
pred_length_sc-2001 -0.03 0.05
pred_length_sc-2002 0.22 0.08
pred_length_sc-2003 -0.06 0.05
pred_length_sc-2004 0.04 0.04
pred_length_sc-2005 0.04 0.05
pred_length_sc-2006 0.31 0.04
pred_length_sc-2007 0.23 0.03
pred_length_sc-2008 0.43 0.03
pred_length_sc-2009 0.42 0.03
pred_length_sc-2010 0.29 0.04
pred_length_sc-2011 0.40 0.15
pred_length_sc-2012 0.50 0.05
pred_length_sc-2013 0.37 0.03
pred_length_sc-2014 0.64 0.13
pred_length_sc-2015 0.65 0.09
pred_length_sc-2016 0.39 0.06
pred_length_sc-2017 0.13 0.05
pred_length_sc-2018 0.48 0.05
pred_length_sc-2019 0.16 0.07
pred_length_sc-2020 0.45 0.05
pred_length_sc-2021 0.02 0.04
pred_length_sc-2022 0.09 0.05
Dispersion parameter: 0.83
Matérn range: 16.79
Spatiotemporal IID SD: 0.54
ML criterion at convergence: -120535.625
See ?tidy.sdmTMB to extract these values as a data frame.## Coef estimate for predator length
pl_M1 <- as.list(M1$sd_report, "Estimate")
pls_M1 <- as.list(M1$sd_report, "Std. Error")
pl_M2 <- as.list(M2$sd_report, "Estimate")
pls_M2 <- as.list(M2$sd_report, "Std. Error")
# the time varying coef:s are in b_rw_t and for M2, both year intercept and pred length effect are in b_rw_t (an array)
ypl_M1 <- data.frame(Year = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year)), coef.est = pl_M1$b_rw_t, coef.se = pls_M1$b_rw_t, model = "M1")
ypl_M2 <- data.frame(Year = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year)), coef.est = pl_M2$b_rw_t[,2,], coef.se = pls_M2$b_rw_t[,2,], model = "M2")
ypl_M12 <- bind_rows(ypl_M1, ypl_M2)
# Conditional effects of length varying with Year
nd_length <- expand.grid(
pred_length_sc = seq(min(df_noz_doy_sel$pred_length_sc), max(df_noz_doy_sel$pred_length_sc), length.out = 50),
Year = unique(df_noz_doy_sel$Year))
nd_length$depth_sc <- 0
nd_length$doy_sc <- 0
p_M1_length <- predict(M1, newdata = nd_length, se_fit = TRUE, re_form = NA)
p_M2_length <- predict(M2, newdata = nd_length, se_fit = TRUE, re_form = NA)
p_M1_length$model = "M1"
p_M2_length$model = "M2"
p_length <- bind_rows(p_M1_length, p_M2_length)
# Coefficient estimate of Year
M1_est <- tidy(M1, effects = "fixed", conf.int = TRUE)$estimate
M1_se <- tidy(M1, effects = "fixed", conf.int = TRUE)$std.error
M2_est <- pl_M2$b_rw_t[,1,]
M2_se <- pls_M2$b_rw_t[,1,]
y_M1 <- data.frame(Year = sort(unique(M1$data$Year)), coef.est = M1_est[1:41], coef.se = M1_se[1:41], model = "M1")
y_M2 <- data.frame(Year = sort(unique(M2$data$Year)), coef.est = M2_est[1:41], coef.se = M2_se[1:41], model = "M2")
# Conditional effect of Year
nd_Year <- data.frame(Year = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year)), depth_sc = 0, pred_length_sc = 0, doy_sc = 0)
p_M1_Year <- predict(M1, newdata = nd_Year, re_form = NA, se_fit = TRUE)
p_M2_Year <- predict(M2, newdata = nd_Year, re_form = NA, se_fit = TRUE)
p_M1_Year$model = "M1"
p_M2_Year$model = "M2"
p_Year <- bind_rows(p_M1_Year, p_M2_Year)
ypl_M12 |>
ggplot(aes(Year, exp(coef.est), color = model)) +
geom_line() +
geom_ribbon(aes(ymin = exp(coef.est-1.96*coef.se), ymax = exp(coef.est+1.96*coef.se)), alpha = 0.3) +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) +
ylab("est for pred.length") 
p_length |>
filter(Year %in% seq(min(Year), max(Year), by = 5)) |>
ggplot(aes(pred_length_sc, exp(est), group = as.factor(Year))) +
facet_wrap(~model) +
geom_line(aes(colour = Year), lwd = 1) #+filter: removed 3,100 rows (79%), 800 rows remaining
geom_ribbon(aes(fill = Year), alpha = 0.1)mapping: fill = ~Year
geom_ribbon: na.rm = FALSE, orientation = NA, outline.type = both
stat_identity: na.rm = FALSE
position_identity bind_rows(y_M1,y_M2) |>
ggplot(aes(Year, exp(coef.est), color = model, linetypre = model)) +
geom_line() +
geom_ribbon(aes(ymin = exp(coef.est-1.96*coef.se), ymax = exp(coef.est+1.96*coef.se)), alpha = 0.3) +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) +
ylab("est for year") # the cf gets very large for the missing Year
p_Year |>
ggplot(aes(Year, exp(est), color = model)) +
geom_line() +
geom_ribbon(aes(ymin = exp(est-1.96*est_se), ymax = exp(est+1.96*est_se)), alpha = 0.3) +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) +
ylim(0,max(exp(p_Year$est) * 1.2)) + # the cf gets very large for the missing Year
ylab("prediction for year")
The difference between M1 and M2 is almost non-existent for the estimates for predator length and small for the predictions (the prediction is smooth for M2 and the cfs are different). The estimate for year is different. What am I missing, how can the predictions be the same but the coeeficient differ so much??
Model alterantives to M1 and M2
M3a1 <- # factor Year * pred_length, skip time varying and missing Years
sdmTMB(
data = df_noz_doy_sel,
formula = fr_tot ~ 0 + factor(Year) * pred_length_sc + s(doy_sc, bs = "cc") + depth_sc,
#time_varying = ~ 1 + pred_length_sc,
#extra_time = missing_Years,
#time = "Year",
mesh = mesh_nozdoy_sel,
spatial = "off",
family = Gamma(link = "log"),
)
M3a2 <- # pred_length * factor Year, skip time varying and missing Years
sdmTMB(
data = df_noz_doy_sel,
formula = fr_tot ~ 0 + pred_length_sc * factor(Year) + s(doy_sc, bs = "cc") + depth_sc,
#time_varying = ~ 1 + pred_length_sc,
#extra_time = missing_Years,
#time = "Year",
mesh = mesh_nozdoy_sel,
spatial = "off",
family = Gamma(link = "log"),
)
M3b <- # smooth on length (k = 3) by Year. Warning: NA/NaN function evaluationWarning: NA/NaN function evaluationWarning: The model may not have converged: non-positive-definite Hessian matrix.
sdmTMB(
data = df_noz_doy_sel,
formula = fr_tot ~ 0 + factor(Year) + s(doy_sc, bs = "cc") + depth_sc + s(pred_length_sc, k = 3, by = Year),
#time_varying = ~ 1,
#extra_time = missing_Years,
#time = "Year",
mesh = mesh_nozdoy_sel,
spatial = "off",
family = Gamma(link = "log"),
)Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient
= tmb_obj$gr, : NA/NaN function evaluation
Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient
= tmb_obj$gr, : NA/NaN function evaluationWarning: The model may not have converged: non-positive-definite Hessian
matrix.M3c <- # mean as random factor instead of fixed
sdmTMB(
data = df_noz_doy_sel,
formula = fr_tot ~ 0 + s(doy_sc, bs = "cc") + depth_sc + (0|Year_f),
time_varying = ~ 0 + pred_length_sc,
extra_time = missing_Years,
time = "Year",
mesh = mesh_nozdoy_sel,
spatial = "off", # It still fits a spatial model
family = Gamma(link = "log"),
)
# The interaction with a smooth effect seemingly does'nt work. This prevents the model from finding Year.
# M3d <- # factor Year * s(pred_length, k = 4).
# sdmTMB(
# data = df_noz_doy_sel_M2,
# formula = fr_tot ~ factor(Year) * s(pred_length_sc, k = 4) + s(doy_sc, bs = "cc") + depth_sc,
# #time_varying = ~ 1 + pred_length_sc,
# #extra_time = missing_Years,
# #time = "Year",
# mesh = mesh_nozdoy_sel_M2,
# spatial = "on",
# family = Gamma(link = "log"),
# )
sanity(M3a1) ✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100sanity(M3b) ✔ Non-linear minimizer suggests successful convergence✖ Non-positive-definite Hessian matrix: model may not have convergedℹ Try simplifying the model, adjusting the mesh, or adding priors✔ No extreme or very small eigenvalues detected✖ `ln_phi` gradient > 0.001ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the modelWarning in sqrt(diag(object$cov.fixed)): NaNs produced✖ `b_j` standard error is NAℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `b_j` standard error is NA
ℹ Try simplifying the model, adjusting the mesh, or adding priors✖ `bs` standard error is NAℹ `bs` is an internal parameter presenting the linear compoment of a smootherℹ It's not uncommon for this parameter to have large standard errors,
which can possibly be ignored if the rest of the model looks OKℹ Try simplifying the model, adjusting the mesh, or adding priorsWarning in sqrt(diag(object$cov.fixed)): NaNs producedWarning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced✔ No standard errors look unreasonably largeWarning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100sanity(M3c)✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeqqnorm(residuals(M3a1))
qqline(residuals(M3a1))
qqnorm(residuals(M3c))
qqline(residuals(M3c))
# Only a and c left for comparison
M3a1Model fit by ML ['sdmTMB']
Formula: fr_tot ~ 0 + factor(Year) * pred_length_sc + s(doy_sc, bs = "cc") +
Formula: depth_sc
Mesh: mesh_nozdoy_sel (isotropic covariance)
Data: df_noz_doy_sel
Family: Gamma(link = 'log')
coef.est coef.se
factor(Year)1982 -3.60 0.03
factor(Year)1983 -3.66 0.02
factor(Year)1984 -3.87 0.03
factor(Year)1985 -3.77 0.04
factor(Year)1986 -3.62 0.03
factor(Year)1987 -3.74 0.04
factor(Year)1988 -3.75 0.04
factor(Year)1989 -3.70 0.05
factor(Year)1990 -3.59 0.04
factor(Year)1991 -3.90 0.06
factor(Year)1993 -3.85 0.07
factor(Year)1994 -3.79 0.08
factor(Year)1995 -4.04 0.06
factor(Year)1996 -4.05 0.07
factor(Year)1997 -4.04 0.11
factor(Year)1998 -4.00 0.11
factor(Year)1999 -3.88 0.04
factor(Year)2000 -3.98 0.04
factor(Year)2001 -4.01 0.06
factor(Year)2002 -3.88 0.11
factor(Year)2003 -3.92 0.05
factor(Year)2004 -4.26 0.05
factor(Year)2005 -3.91 0.05
factor(Year)2006 -4.26 0.03
factor(Year)2007 -4.01 0.03
factor(Year)2008 -4.29 0.03
factor(Year)2009 -4.21 0.03
factor(Year)2010 -4.22 0.06
factor(Year)2012 -4.35 0.04
factor(Year)2013 -4.36 0.03
factor(Year)2014 -4.32 0.14
factor(Year)2015 -4.22 0.08
factor(Year)2016 -4.27 0.06
factor(Year)2017 -4.33 0.05
factor(Year)2018 -4.15 0.04
factor(Year)2019 -4.55 0.08
factor(Year)2020 -4.25 0.05
factor(Year)2021 -4.23 0.03
factor(Year)2022 -4.43 0.06
pred_length_sc 0.40 0.03
depth_sc 0.00 0.01
factor(Year)1983:pred_length_sc -0.05 0.04
factor(Year)1984:pred_length_sc -0.13 0.04
factor(Year)1985:pred_length_sc -0.18 0.05
factor(Year)1986:pred_length_sc -0.26 0.05
factor(Year)1987:pred_length_sc -0.28 0.05
factor(Year)1988:pred_length_sc -0.21 0.06
factor(Year)1989:pred_length_sc -0.38 0.06
factor(Year)1990:pred_length_sc -0.27 0.05
factor(Year)1991:pred_length_sc -0.09 0.06
factor(Year)1993:pred_length_sc -0.53 0.07
factor(Year)1994:pred_length_sc -0.35 0.09
factor(Year)1995:pred_length_sc -0.43 0.06
factor(Year)1996:pred_length_sc -0.62 0.07
factor(Year)1997:pred_length_sc -0.03 0.09
factor(Year)1998:pred_length_sc -0.27 0.08
factor(Year)1999:pred_length_sc -0.38 0.05
factor(Year)2000:pred_length_sc -0.33 0.06
factor(Year)2001:pred_length_sc -0.27 0.06
factor(Year)2002:pred_length_sc -0.13 0.09
factor(Year)2003:pred_length_sc -0.32 0.05
factor(Year)2004:pred_length_sc -0.35 0.05
factor(Year)2005:pred_length_sc -0.27 0.05
factor(Year)2006:pred_length_sc 0.04 0.05
factor(Year)2007:pred_length_sc -0.21 0.04
factor(Year)2008:pred_length_sc 0.04 0.04
factor(Year)2009:pred_length_sc -0.09 0.04
factor(Year)2010:pred_length_sc -0.11 0.05
factor(Year)2012:pred_length_sc 0.08 0.05
factor(Year)2013:pred_length_sc -0.07 0.04
factor(Year)2014:pred_length_sc 0.44 0.27
factor(Year)2015:pred_length_sc 0.09 0.09
factor(Year)2016:pred_length_sc -0.06 0.07
factor(Year)2017:pred_length_sc -0.28 0.05
factor(Year)2018:pred_length_sc 0.13 0.06
factor(Year)2019:pred_length_sc -0.33 0.08
factor(Year)2020:pred_length_sc 0.01 0.06
factor(Year)2021:pred_length_sc -0.27 0.05
factor(Year)2022:pred_length_sc -0.29 0.06
Smooth terms:
Std. Dev.
sds(doy_sc) 0.1
Dispersion parameter: 0.77
ML criterion at convergence: -119318.108
See ?tidy.sdmTMB to extract these values as a data frame.M3b # Warning: Smoother fixed effect matrix names could not be retrieved. Thi smakes the smoother effect Model fit by ML ['sdmTMB']
Formula: fr_tot ~ 0 + factor(Year) + s(doy_sc, bs = "cc") + depth_sc +
Formula: s(pred_length_sc, k = 3, by = Year)
Mesh: mesh_nozdoy_sel (isotropic covariance)
Data: df_noz_doy_sel
Family: Gamma(link = 'log')Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs producedWarning: Smoother fixed effect matrix names could not be retrieved.
This could be from setting m != 2 in s() or other unanticipated s() arguments.
This does not affect model fitting.
We'll use generic covariate names ('scovariate') here intead.Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
coef.est coef.se
factor(Year)1982 -3.60 NaN
factor(Year)1983 -3.64 NaN
factor(Year)1984 -3.85 NaN
factor(Year)1985 -3.78 NaN
factor(Year)1986 -3.66 NaN
factor(Year)1987 -3.78 NaN
factor(Year)1988 -3.77 NaN
factor(Year)1989 -3.86 NaN
factor(Year)1990 -3.64 NaN
factor(Year)1991 -3.85 NaN
factor(Year)1993 -3.96 NaN
factor(Year)1994 -3.80 NaN
factor(Year)1995 -4.08 NaN
factor(Year)1996 -4.24 NaN
factor(Year)1997 -3.86 NaN
factor(Year)1998 -3.83 NaN
factor(Year)1999 -3.82 NaN
factor(Year)2000 -3.96 NaN
factor(Year)2001 -3.86 NaN
factor(Year)2002 -3.86 NaN
factor(Year)2003 -3.76 NaN
factor(Year)2004 -4.04 NaN
factor(Year)2005 -3.83 NaN
factor(Year)2006 -4.27 NaN
factor(Year)2007 -3.99 NaN
factor(Year)2008 -4.28 NaN
factor(Year)2009 -4.15 NaN
factor(Year)2010 -4.12 NaN
factor(Year)2012 -4.36 NaN
factor(Year)2013 -4.36 NaN
factor(Year)2014 -4.55 NaN
factor(Year)2015 -4.30 NaN
factor(Year)2016 -4.25 NaN
factor(Year)2017 -4.25 NaN
factor(Year)2018 -4.20 NaN
factor(Year)2019 -4.36 NaN
factor(Year)2020 -4.31 NaN
factor(Year)2021 -4.19 NaN
factor(Year)2022 -4.29 NaN
depth_sc -0.01 0.01
scovariate-1 0.00 0.00
scovariate-2 0.00 NaN
Smooth terms:
Std. Dev.
sds(doy_sc) 0.11
sds(pred_length_sc):Year) 0.00
Dispersion parameter: 0.77
ML criterion at convergence: -119107.484
See ?tidy.sdmTMB to extract these values as a data frame.
**Possible issues detected! Check output of sanity().**M3cSpatiotemporal model fit by ML ['sdmTMB']
Formula: fr_tot ~ 0 + s(doy_sc, bs = "cc") + depth_sc + (0 | Year_f)
Mesh: mesh_nozdoy_sel (isotropic covariance)
Time column: Year
Data: df_noz_doy_sel
Family: Gamma(link = 'log')
coef.est coef.se
depth_sc 0 0.01
Smooth terms:
Std. Dev.
sds(doy_sc) 0.06
Random intercepts:
Std. Dev.
Year_f 4.15
Time-varying parameters:
coef.est coef.se
pred_length_sc-1982 0.36 0.04
pred_length_sc-1983 0.27 0.03
pred_length_sc-1984 0.23 0.04
pred_length_sc-1985 0.09 0.04
pred_length_sc-1986 0.10 0.04
pred_length_sc-1987 0.06 0.04
pred_length_sc-1988 0.17 0.05
pred_length_sc-1989 0.05 0.05
pred_length_sc-1990 0.10 0.04
pred_length_sc-1991 0.26 0.05
pred_length_sc-1992 0.09 0.15
pred_length_sc-1993 -0.12 0.07
pred_length_sc-1994 -0.04 0.08
pred_length_sc-1995 0.00 0.05
pred_length_sc-1996 -0.12 0.07
pred_length_sc-1997 0.26 0.08
pred_length_sc-1998 0.10 0.09
pred_length_sc-1999 -0.01 0.04
pred_length_sc-2000 -0.12 0.06
pred_length_sc-2001 -0.03 0.05
pred_length_sc-2002 0.21 0.08
pred_length_sc-2003 -0.06 0.06
pred_length_sc-2004 0.02 0.04
pred_length_sc-2005 0.06 0.05
pred_length_sc-2006 0.31 0.04
pred_length_sc-2007 0.23 0.03
pred_length_sc-2008 0.43 0.03
pred_length_sc-2009 0.42 0.03
pred_length_sc-2010 0.28 0.04
pred_length_sc-2011 0.39 0.15
pred_length_sc-2012 0.49 0.05
pred_length_sc-2013 0.37 0.03
pred_length_sc-2014 0.64 0.14
pred_length_sc-2015 0.67 0.09
pred_length_sc-2016 0.39 0.06
pred_length_sc-2017 0.12 0.05
pred_length_sc-2018 0.48 0.05
pred_length_sc-2019 0.14 0.07
pred_length_sc-2020 0.46 0.05
pred_length_sc-2021 0.02 0.04
pred_length_sc-2022 0.08 0.06
Dispersion parameter: 0.83
Matérn range: 17.11
Spatiotemporal IID SD: 0.54
ML criterion at convergence: -120413.922
See ?tidy.sdmTMB to extract these values as a data frame.# estimates for Year effects
M3a1_est <- tidy(M3a1, effects = "fixed", conf.int = TRUE)$estimate
M3a1_se <- tidy(M3a1, effects = "fixed", conf.int = TRUE)$std.error
y_M3a <- data.frame(Year = sort(unique(M3a1$data$Year)), coef.est = M3a1_est[1:39], coef.se = M3a1_se[1:39], model = "M3a")
y_M3c <- data.frame(Year = sort(unique(df_noz_doy_sel$Year)), coef.est = tidy(M3c, effects = "ran_vals", conf.int = TRUE)$estimate, coef.se = tidy(M3c, effects = "ran_vals", conf.int = TRUE)$std.error, model = "M3c")
y_M3a$model = "M3a"
y_M3c$model = "M3c"
y_M3 <- bind_rows(y_M3a,y_M3c)
# Conditional effect of Year. Cant predict over Year when it is a random effect somehow
# nd_M3a_Year <- data.frame(Year = sort(unique(M3a$data$Year)), depth_sc = 0, pred_length_sc = 0, doy_sc = 0)
# nd_M3c_Year <- data.frame(Year = sort(unique(M3c$data$Year)), depth_sc = 0, pred_length_sc = 0, doy_sc = 0)
# p_M3a_Year <- predict(M3a, newdata = nd_M3a_Year, re_form = NA, se_fit = TRUE)
# p_M3c_Year <- predict(M3c, newdata = nd_M3c_Year, re_form = NA, se_fit = TRUE) # cant predict over Year...
# Predator length effects varying with Year
# M3a
ypl_M3a1 <- data.frame(Year = sort(unique(M3a1$data$Year)), coef.est = c(M3a1_est[1] + M3a1_est[40], M3a1_est[1] + M3a1_est[40] + M3a1_est[42:79]), coef.se = c(M3a1_se[1], M3a1_se[42:79]), model = "M3a1")
M3a2_est <- tidy(M3a2, effects = "fixed", conf.int = TRUE)$estimate
ypl_M3a2 <- data.frame(Year = sort(unique(M3a2$data$Year)), coef.est = c(M3a2_est[1] + M3a2_est[2], M3a2_est[1] + M3a2_est[2] + M3a2_est[42:79]), model = "M3a2")
ypl_comp_M3a <- bind_rows(ypl_M3a1[,c(1,2,4)], ypl_M3a2)
# M3c
pl_M3c <- as.list(M3c$sd_report, "Estimate")
pls_M3c <- as.list(M3c$sd_report, "Std. Error")
M3c_est <- tidy(M3c, effects = "ran_vals", conf.int = TRUE)$estimate # random Year effect
M3c_se <- tidy(M3c, effects = "ran_vals", conf.int = TRUE)$std.error
ypl_M3c <- data.frame(Year = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year)), coef.est = M3c_est[1] + pl_M3c$b_rw_t, coef.se = pls_M3c$b_rw_t, model = "M3c")
# plot comparison of M3a1 and 2 of the effect of pred length varying with year
ypl_comp_M3a |>
ggplot(aes(Year, exp(coef.est), color = model, linetype = model)) +
geom_line() +
ggtitle(" order of interactive terms (a*b vs b*a) makes no difference")
# plot Year coef for a and c
y_M3 |>
ggplot(aes(Year, exp(coef.est), color = model)) +
geom_line() +
geom_ribbon(aes(ymin = exp(coef.est-1.96*coef.se), ymax = exp(coef.est+1.96*coef.se)), alpha = 0.3) +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 20)) +
ylab("random estimate of Year")
# compare M3a and c of the effect if pred length varying with year
bind_rows(ypl_M3c, ypl_M3a1) |>
ggplot(aes(Year, exp(coef.est), color = model)) +
geom_line() +
geom_ribbon(aes(ymin = exp(coef.est-1.96*coef.se), ymax = exp(coef.est+1.96*coef.se)), alpha = 0.3) +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) +
ylab("est for pred.length") 
In the plots above, M3a has the intercept for the first year (1982) and the fixed effect of pred_length_sc added to the interaction effect of year and pred_length_sc. M3c has the random effect of year added to the effect of predator length varying with year but no overall effect of predator length. This is a bit odd as the predator length effect is large (est 0.4) in M3a. What am I missing..
# Compare M3 with M1 and 2
bind_rows(ypl_M3a1,ypl_M3c) |>
bind_rows(ypl_M12) |>
ggplot(aes(Year, exp(coef.est), color = model, linetype = model)) +
geom_line() +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) +
ylab("est for pred.length")
The intercepts may be different…
# M1
# formula = fr_tot ~ 0 + as.factor(Year) + s(doy_sc, bs = "cc") + depth_sc,
# time_varying = ~ 0 + pred_length_sc, # implies time varying random walk intercept
# M2
# formula = fr_tot ~ 0 + s(doy_sc, bs = "cc") + depth_sc,
# time_varying = ~ 1 + pred_length_sc,
# M3a
# formula = fr_tot ~ 0 + factor(Year) * pred_length_sc + s(doy_sc, bs = "cc") + depth_sc,
# M3c
# formula = fr_tot ~ 0 + s(doy_sc, bs = "cc") + depth_sc + (0|Year_f),
# time_varying = ~ 0 + pred_length_sc,
# I remove the year intercepts, i.e. 1982 for M3a (but keep the predator length effect) and the random year effect for M3c..
# ypl_M3a1 <- data.frame(Year = sort(unique(M3a1$data$Year)), coef.est = c(M3a1_est[1] + M3a1_est[40], M3a1_est[1] + M3a1_est[40] + M3a1_est[42:79]), coef.se = c(M3a1_se[1], M3a1_se[42:79]), model = "M3a1")
ypl_M3aa <- data.frame(Year = sort(unique(M3a1$data$Year)), coef.est = c(M3a1_est[40], M3a1_est[40] + M3a1_est[42:79]), coef.se = c(M3a1_se[1], M3a1_se[42:79]), model = "M3a1")
# ypl_M3c <- data.frame(Year = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year)), coef.est = M3c_est[1] + pl_M3c$b_rw_t, coef.se = pls_M3c$b_rw_t, model = "M3c")
ypl_M3cc <- data.frame(Year = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year)), coef.est = pl_M3c$b_rw_t, coef.se = pls_M3c$b_rw_t, model = "M3c")
bind_rows(ypl_M3aa,ypl_M3cc) |>
bind_rows(ypl_M12) |>
ggplot(aes(Year, exp(coef.est), color = model, linetype = model)) +
geom_line() +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) +
ylab("est for pred.length") 
By removing the year intercepts, i.e. 1982 for M3a and the random year effect for M3c, we get the coeffcients for the effects of predator length varying with year on the same scale as M1 and M2.
Summary
- M3a - Specifying pred_length_sc\(*\)factor(Year) or factor(Year)\(*\)pred_length_sc doesn’t matter. A fixed effect interactiopn makes for somewhat different coefficient estimates for predator length over year compared to having this effect in time_varying.
- M3b - The effect of a smoother on predator length varying with Year (i.e. by = Year) cannot be estimated by the model( ” Smoother fixed effect matrix names could not be retrieved”). The smoother term can’t be put on time_varying (the function doesn’t recognize the s() function). There are more issue with the model (see above) and we need to look into this if we want to do this model.I’ve dropped this option.
- M3c - year as a randoom effect instead of a time_varying intercept makes small differences.
- M3d - An interaction between Year and a smooth effect on pred_length_sc seemingly doesn’t work. This prevents sdmTMB from finding Year. Therefore I’ve dropped this option.
Independent of model structure, the models mantain the same estimate for the predator length effects varying over Year. The Year effect may differ however
I also NOTEd that, if time_varying is on but spatial is off, it fits a spatial model and asks for xy coordinates when I use predict().
Appendix 1. Variation and time_varying intercepts vs fixed factor year
# Pred length i fixed effect. test på 10 år. Om m2 smoothen är, t_v contra fator år har lite skilllnad men här har vi. Intercpetrn er . testa genom längs om en fixe deefkt om tv fixed ekfkket kan fånaga annual fluctuatio
# M1 varies but not M2
Mtesta <-
sdmTMB(
data = df_noz_doy_sel |> filter(Year_f %in% 2012:2022),
formula = fr_tot ~ 0 + s(doy_sc, bs = "cc") + depth_sc,
time_varying = ~ 1 + pred_length_sc,
time = "Year", # for spatiotemporal and time_varying
#mesh = mesh_nozdoy_sel,
spatiotemporal = "off",
spatial = "off",
family = Gamma(link = "log"),
)filter: removed 29,705 rows (75%), 9,649 rows remainingMtestb <-
sdmTMB(
data = df_noz_doy_sel |> filter(Year_f %in% 2012:2022),
formula = fr_tot ~ 0 + Year_f + s(doy_sc, bs = "cc") + depth_sc,
time_varying = ~ 0 + pred_length_sc,
#extra_time = missing_Years, # to fill in empty Years
time = "Year", # for spatiotemporal and time_varying
#mesh = mesh_nozdoy_sel,
spatial = "off",
spatiotemporal = "off",
family = Gamma(link = "log"),
)filter: removed 29,705 rows (75%), 9,649 rows remainingMtestc <-
sdmTMB(
data = df_noz_doy_sel |> filter(Year_f %in% 2012:2022),
formula = fr_tot ~ 0 + pred_length_sc + s(doy_sc, bs = "cc") + depth_sc,
time_varying = ~ 1,
#extra_time = missing_Years, # to fill in empty Years
time = "Year", # for spatiotemporal and time_varying
#mesh = mesh_nozdoy_sel,
spatial = "off",
spatiotemporal = "off",
family = Gamma(link = "log"),
)filter: removed 29,705 rows (75%), 9,649 rows remainingMtestd <-
sdmTMB(
data = df_noz_doy_sel |> filter(Year_f %in% 2012:2022),
formula = fr_tot ~ 0 + Year_f + pred_length_sc + s(doy_sc, bs = "cc") + depth_sc,
time_varying = ~ 0,
time = "Year", # for spatiotemporal and time_varying
#mesh = mesh_nozdoy_sel,
spatiotemporal = "off",
spatial = "off",
family = Gamma(link = "log"),
)filter: removed 29,705 rows (75%), 9,649 rows remainingnd_Year <- data.frame(Year = Mtesta$data$Year, depth_sc = 0, pred_length_sc = 0, doy_sc = 0)
p_Mtesta_Year <- predict(Mtesta, newdata = nd_Year, re_form = NA, se_fit = TRUE)
nd_Year <- data.frame(Year = Mtestb$data$Year, Year_f = Mtestb$data$Year_f, depth_sc = 0, pred_length_sc = 0, doy_sc = 0)
p_Mtestb_Year <- predict(Mtestb, newdata = nd_Year, re_form = NA, se_fit = TRUE)
nd_Year <- data.frame(Year = Mtestc$data$Year, depth_sc = 0, pred_length_sc = 0, doy_sc = 0)
p_Mtestc_Year <- predict(Mtestc, newdata = nd_Year, re_form = NA, se_fit = TRUE)
nd_Year <- data.frame(Year = Mtestd$data$Year, Year_f = Mtestd$data$Year_f, depth_sc = 0, pred_length_sc = 0, doy_sc = 0)
p_Mtestd_Year <- predict(Mtestd, newdata = nd_Year, re_form = NA, se_fit = TRUE)
p_Mtesta_Year$model = "a"
p_Mtestb_Year$model = "b"
p_Mtestc_Year$model = "c"
p_Mtestd_Year$model = "d"
p_Year <- bind_rows(p_Mtesta_Year, p_Mtestb_Year,p_Mtestc_Year, p_Mtestd_Year)
p_Year |>
ggplot(aes(Year, exp(est), color = model)) +
geom_line() +
geom_ribbon(aes(ymin = exp(est-1.96*est_se), ymax = exp(est+1.96*est_se)), alpha = 0.1) +
scale_x_continuous(breaks = seq(min(df_noz_doy_sel$Year), max(df_noz_doy_sel$Year), by = 10)) +
ylab("prediction for year")
Appendix 2. Model with a smoother on length
This causes issues that I so far havent figured out how to solve
M0 <-
sdmTMB(
data = df_noz_doy_sel |> filter(Year > 2012),
formula = fr_tot ~ 0 + Year_f + s(pred_length_sc, m = 2, by = Year_f) + depth_sc,
spatiotemporal = "off",
spatial = "off",
family = Gamma(link = "log"),
)filter: removed 30,660 rows (78%), 8,694 rows remainingsanity(M0)✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✖ `ln_smooth_sigma` standard error may be largeℹ Try simplifying the model, adjusting the mesh, or adding priors✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100M0Model fit by ML ['sdmTMB']
Formula: fr_tot ~ 0 + Year_f + s(pred_length_sc, m = 2, by = Year_f) +
Formula: depth_sc
Mesh: NULL (isotropic covariance)
Data: filter(df_noz_doy_sel, Year > 2012)
Family: Gamma(link = 'log')
coef.est coef.se
Year_f2013 -4.51 0.03
Year_f2014 -4.55 0.10
Year_f2015 -4.41 0.07
Year_f2016 -4.40 0.06
Year_f2017 -4.49 0.05
Year_f2018 -4.32 0.05
Year_f2019 -4.70 0.09
Year_f2020 -4.33 0.06
Year_f2021 -4.35 0.03
Year_f2022 -4.54 0.08
depth_sc 0.03 0.01
spred_length_sc):Year_f2013 0.78 0.86
spred_length_sc):Year_f2014 -0.71 0.25
spred_length_sc):Year_f2015 -0.41 0.08
spred_length_sc):Year_f2016 -0.28 0.05
spred_length_sc):Year_f2017 1.16 1.52
spred_length_sc):Year_f2018 1.09 2.19
spred_length_sc):Year_f2019 1.10 2.01
spred_length_sc):Year_f2020 0.10 1.06
spred_length_sc):Year_f2021 0.81 0.51
spred_length_sc):Year_f2022 1.10 1.59
Smooth terms:
Std. Dev.
sds(pred_length_sc):Year_f2013) 6.54
sds(pred_length_sc):Year_f2014) 0.00
sds(pred_length_sc):Year_f2015) 0.00
sds(pred_length_sc):Year_f2016) 0.00
sds(pred_length_sc):Year_f2017) 9.50
sds(pred_length_sc):Year_f2018) 13.57
sds(pred_length_sc):Year_f2019) 12.64
sds(pred_length_sc):Year_f2020) 6.35
sds(pred_length_sc):Year_f2021) 17.43
sds(pred_length_sc):Year_f2022) 9.23
Dispersion parameter: 0.65
ML criterion at convergence: -30622.212
See ?tidy.sdmTMB to extract these values as a data frame.
**Possible issues detected! Check output of sanity().**nd_length <- expand.grid(
pred_length_sc = seq(min(df_noz_doy_sel$pred_length_sc), max(df_noz_doy_sel$pred_length_sc), length.out = 50),
Year_f = factor(2013:2022))
nd_length$depth_sc <- 0
nd_length$doy_sc <- 0
p_M0_length <- predict(M0, newdata = nd_length, se_fit = TRUE, re_form = NA)
p_M0_length |>
ggplot(aes(pred_length_sc, exp(est), color = Year_f)) +
geom_line() #+
geom_ribbon(aes(fill = Year), alpha = 0.1)mapping: fill = ~Year
geom_ribbon: na.rm = FALSE, orientation = NA, outline.type = both
stat_identity: na.rm = FALSE
position_identity